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The 2D time dependent solution of thin accretion disk in a close binary system have 
been presented on the equatorial plane around the Schwarzschild black hole. To do 
that, the special part of the General Relativistic Hydrodynamical(GRH) equations 
are solved using High Resolution Shock Capturing (HRSC) schemes. The spiral 
shock waves on the accretion disk are modeled using perfect fluid equation of state 
with adiabatic indices 7 = 1.05,1.2 and 5/3. The results show that the spiral shock 
waves are created for gammas except the case 7 = 5/3. These results consistent with 
results from Newtonian hydrodynamic code except close to black hole. Newtonian 
approximation does not give good solution while matter closes to black hole. Our 
simulations illustrate that the spiral shock waves are created close to black hole and 
the location of inner radius of spiral shock wave is around 10M and it depends on 
the specific heat rates. We also find that the smaller 7 is the more tightly the spiral 
winds. 

Keywords: General Relativity, Hydrodynamics, Numerical Relativity Black Hole, Accretion 
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I. INTRODUCTION 

Rotating accretion disk around compact objects is an important problem in astrophysics, 
such as neutron stars and black holes which involve mass transfer from one object to another. 
Shock waves in a rotating accretion disk onto compact objects transfer the gravitation energy 
to the radiation energy which is observer by different X-ray observatory satellite, such as 
Chandra. Fluctuations and oscillating of radiation emitted from these systems constantly 
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remind us of time variations of the dynamical quantities. These variation can happen on 
times scale of a few microseconds to few years. In order to understand these kinds of events 
we have started doing some numerical simulation and looking for shock waves on an accretion 
disk. 

In binary systems in which a compact primary star accretes, through an accretion disk 
from a lobe-filling secondary star, the tidal interaction with the companion can results in 
the formation of a two-armed spiral structure(Dgani et al. 1994, Godon et al. 1997). These 
spiral waves can transport angular momentum in the accretion disk which has been suggested 
by analogy to the galactic dynamics context (Rozyczka et al. 1993). 

The problem of accretion disk on to black hole has been previously analyzed numerically 
by Hawley et al. (1985) in general relativity. They have also suggested shock waves on the 
rotating accretion disk. After that so many astrophysicist worked on this kinds of problem 
using with Newtonian (pseudo-Newtonian gravitational potential approximation is used) and 
general relativistic hydrodynamics. In Chakabarti et al. (1993), They have compared the 
analytic and numerical studies of shock wave on inviscid accretion disk flows using with the 
smoothed particle hydrodynamics (SPH) code. They have believed that shocks waves could 
be common in accretion disk and standing shock waves can only be produced in accretion disk 
around the black hole. Near to black hole density of accretion disk in the sub-Keplerian flows 
is higher than density of accretion disk in Keplerian flows the consequence of the presence 
of centrifugal barrier, which is smaller at sub-Keplerian flow. Because of this high density 
standing shock can form in accretion disk. This high density flow intercepts soft photons 
from a cold Keplerian disk and reprocesses them to form high energy X-rays (Chakabarti et 
al. 1997). In Lanzafame et al. (1998), they have used SPH code with viscosity to look at 
the shock waves on an accretion disk at parameter range. They have found that if viscosity 
parameter is less than critical values, shock can form. If it is bigger than critical values, shock 
wave disappears. In the intermediate viscosity, the disk oscillates in the viscous time scales. 
Nonaxisymmetric shock waves are found on the equatorial plane by Molteni et al. (1999). 
They have used SPH and Eulerian type with TVD codes to look for nonaxisymmetric shock 
wave applying a perturbation on an accretion disk at different parameter range, which are 
specific angular momentum and internal energy. They have concluded that shock waves are 
found by Chakabarti with perturbation in a rotating inviscid accretion flow are generally 
unstable to azimuthal perturbation, but the instability are taken care of at low level and 
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new stable asymmetric accretion disk is developed with a strong shock rotating steadily. 

In Makita et al. (2000), first, they have reviewed the spiral shock wave in 2 D and 3D 
and showed their results to look for consistency with literature. One of the main problem 
in accretion disk is the mechanism of angular momentum transport. One of the reason for 
angular momentum transport is a-disk model. The viscosity is supposed to transform the 
angular momentum and it can be produce shock waves depends on the that parameter. The 
another way of the transporting angular momentum is the spiral shock waves. The first 
convincing evidence of spiral shock wave in a accretion disk is observed by Steeghs et al. 
(1997). They used the technique which is known as Doppler temography to observe spiral 
structure in the accretion disk of the eclipsing dwarf nova binary IP Peg at the outburst 
phase. 

Here we do perturbation onto satirically symmetric steady state accretion disk around 
the black hole. Most of the numerical calculation for accretion disk are done by Newtonian 
hydrodynamical code using relativistic approximation on it. But our code fully inviscid 
general relativistic hydrodynamical which is used High Resolution Shock Capturing Scheme 
(HRSC). General relativistic code gives us more detail explanations when the fluid flow 
closes to black hole. 


II. FORMULATION 

The GRH equations in references Font et al. (2000) and Donat et al. (1998), written 
in the standard covariant form, consist of the local conservation laws of the stress-energy 
tensor T^ v and the matter current density JG 


V m T^ = 0, Vm^ = 0- (!) 

Greek indices run from 0 to 3, Latin indices from 1 to 3, and units in which the speed of 
light c = 1 are used. 

Defining the characteristic waves of the general relativistic hydrodynamical equations 
is not trivial with imperfect fluid stress-energy tensor. We neglect the viscosity and heat 
conduction effects. This defines the perfect fluid stress-energy tensor. We use this stress- 
energy tensor to derive the hydrodynamical equations. With this perfect fluid stress-energy 
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tensor, we can solve some problems which are solved by the Newtonian hydrodynamics 
with viscosity, such as those involving angular momentum transport and shock waves on an 
accretion disk, etc. Entropy for perfect fluid is conserved along the fluid lines. The stress 
energy tensor for a perfect fluid is given as 

T> 1U = phuW + Pg^. (2) 

A perfect fluid is a fluid that moves through spacetime with a 4-velocity tA which may vary 
from event to event. It exhibits a density of mass p and isotropic pressure P in the rest 
frame of each fluid element, h is the specific enthalpy, defined as 

h — 1 + e H-. (3) 

P 

Here e is the specific internal energy. The equation of state might have the functional form 
P = P(p, e). The perfect gas equation of state, 

p = (r - 1 )pe, (4) 

is such a functional form. 

The conservation laws in the form given in Eq. © are not suitable for the use in advanced 
numerical schemes. In order to carry out numerical hydrodynamic evolutions such as those 
reported in Font et al. (2000), and to use HRSC methods, the hydrodynamic equations 
after the 3+1 split must be written as a hyperbolic system of first order flux conservative 
equations. We write Eq. (JTJ) in terms of coordinate derivatives, using the coordinates (x° = 
t, x 1 , x 2 , x 3 ). Eq. © is projected onto the basis {nfy (^fy), where rp is a unit timclike 
vector normal to a given hypersurface. After a straightforward calculation, we get (see Font 
et al. 2000), 


d t U + diF* — S, (5) 

where dt = d/dt and di = d/dx 1 . This basic step serves to identify the set of unknowns, the 
vector of conserved quantities U, and their corresponding fluxes F{U). With the equations in 
conservation form, almost every high resolution method devised to solve hyperbolic systems 
of conservation laws can be extended to GRH. 
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The evolved state vector U consists of the conservative variables ( D,Sj,r ) which are 
conserved variables for density, momentum and energy respectively; in terms of the primitive 
variables (p,v\e), this becomes (Font et al. 2000) 



( D N 


1 ^/Wp N 

u = 

s* 

= 

y/^phW 2 Vj 


\ T ) 


y ^7 (phW 2 — P — Wp) j 


( 6 ) 


Here 7 is the determinant of the 3-metric 7 ^, Vj is the fluid 3-velocity, and W is the Lorentz 
factor, 


W = au° = (1 — 7 ijV l v j ) 1 
The flux vectors F l are given by [? ] 


(7) 


/ 


F l = 


a {v* - 


\ 


( 8 ) 


<*{(?*-fflSj + yfyPfy 

\ a {( yi - ) r + V^ yip } ) 

The spatial components of the 4-velocity u l are related to the 3-velocity by the following 
formula: u l = W{v l — j3 l /a). a and f3 l are the lapse function and the shift vector of the 
spacetime respectively. The source vector S is given by Font et al. (2000) 


( 0 1 

s = a^T^g m 

y OLy/liT^dpOL - olT^TIJ) J 
where T^ v is the 4-dimensional Christoffel symbol 


(9) 


l'"„ = + a„9«J - (10) 

The numerical solution of general relativistic hydrodynamical equations and technique 
used are explained in detail our first paper Donmez (2003) which gives great detail about 
formulations, numerical scheme, technique, numerical solution of GRH equations, Adaptive- 
Mesh Refinement (AMR) and solution of special relativistic test problem. 
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III. GENERAL RELATIVISTIC HYDRODYNAMICAL TEST PROBLEM 

In this section, the Schwarzschild geometry is introduced in spherical coordinates to 
define sources for the general relativistic hydrodynamical equations. Then the accretion 
disk problems are numerically modeled , which have been analytically analyzed, to test the 
full GRH code in 2D in the equatorial plane. 


A. Schwarzschild Black Hole 


The Schwarzschild solution determined by the mass M gives the geometry in outside of 
a spherical star or black hole. The Schwarzschild spacetime metric in spherical coordinates 

is 



r 


r 


r 2 dd 2 + r 2 sin 2 6d(f) 2 . 


( 11 ) 


It behaves badly near r = 2M; there the first term becomes zero and the second term 


becomes infinite in Eq. HU). That radius r = 2 M is called the Schwarzschild radius or the 
Schwarzschild horizon. 


The spacetime metric for this line element is 


9 nv 


/-(1-2M) 0 0 0 \ 

0 ( l -^)" 1 0 0 

0 0 r 2 0 


( 12 ) 


V o 


0 0 r 2 sin 2 9 ) 


The lapse function and shift vector for this metric is 


I3 r = 0.0, !i" = 0.0, a" = 0.0, 



( 13 ) 


r 
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B. The Source Terms For Schwarzschild Coordinates 


The gravitational sources for the GRH equations are given by Eq. ©. In order to compute 
the sources in Schwarzschild coordinates for different conserved variables, Eq. © can be 
rewritten as, 


/ 


\ 


0 


a^ T ^g vc Ylr 

a^ T ^g va Ylo 

V - aT^Y^jJ) ) 


(14) 


It is seen in Eq.© that the source for conserved density, D, is zero but the other sources 
depend on the components of the stress energy tensor, Christoffel symbols, and 4-metric. 
After doing some straightforward calculations, the sources can be rewritten in Schwarzschild 
coordinates for each conserved variable with the following form: 

The source for the momentum equation in the radial direction is 


oi\f^T iLV 'gvaY^j. — -jOL\f^j{T tt d r gtt + T rr d r g rr + 

T oe d r goo + T^d r g u ). 


(15) 


The source for the momentum equation in the 9 direction is 


r;» = \aVi T **deau- 


(16) 


The source for the momentum equation in the (j) direction is 


a \/lg ua Y — 0.0. 


(17) 


The source for the energy equation is 


a^T^a - aT<“ ry = 
ay/g(T H d r a - aT rt g tt d r g tt ). 


( 18 ) 




The non-zero components of the stress-energy tensor in Schwarzschild coordinates can be 
computed by Eq. ©; they are 


,, W 2 
T a = ph— + Pg tt 
cr 

T rr = phW 2 {v r ) 2 + Pg rr 
T ee = phW 2 (v e f + Pg ee 
T H = phW 2 {y*) 2 + Pg H 


(19) 


C. Geodesics Flows 

As a general relativistic test problem, the accretion of dust particles onto a black hole 
are solved. The exact solution for pressureless dust is given in Appendix El This problem 
is numerically analyzed in 2D in spherical coordinates at constant 6 = 7t/2, which is the 
equatorial plane, so the spatial numerical domain is the (r, 0) plane. In this calculation,2.4 < 
r < 20 and 0 < 0 < 27T are used for the computational domain. The initial conditions for all 
variables are chosen to have negligible values except the outer boundary (r = 20M) where 
gas is continuously injected radially with the analytic density and velocity. Throughout the 
calculation, whenever values at the outer boundary are needed the analytic values are used. 
The code is tun until a steady state solution is reached, using outflow boundary conditions at 
r = 2.4, inflow boundary conditions at r = 20 and the periodic boundary in the 0 direction. 
It is found that the resulting numerical solution does not develop any angular dependence 
during the simulation. 

In Fig. [U the rest-mass density p, absolute velocity v = {i^Vi) 1 ^ 2 and radial velocity v r 
are plotted as a function of radial coordinate at a fixed angular position. The numerical 
solution agrees well with the analytic solution. The convergence test are also made on this 
problem to test the behavior of the source terms in the GR Hydro code. Some convergence 
tests with the SR Hydro code are conducted to confirm the second order convergent (Donmez 
2004). In this case, we are looking for convergence rate with source terms. The analytic 
values of the accretion problem are used as initial conditions, and the computational domain 
is chosen from r min = 10M to r max = 30M. The convergence results are given at Tabled It 
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TABLE I: L\ norm error and convergence factors are given for different resolutions. 


Convergence Test 

# of points 

# of time step 

L i norm error 

Convergence factor 

32 

1 

1.675 E - 5 


64 

2 

4.0811A -6 

4.105 

128 

4 

9.6456A - 7 

4.23 

256 

8 

2.1372A - 7 

4.51 


is noticed that code gives roughly second order convergence. 

The conservation form of general relativistic hydro dynamical equations are solved. It is 
expected that conserved variables must be conserved to machine accuracy, r-u 10“ 16 . The 
checking the results of conservation variables in the numerical test problem shows that these 
variables conserved to machine accuracy. 

In this part of same test problem we do an AMR test. For uniform grid runs, the amount 
of time it takes to reach a steady state solution increases with resolution. We carried out 
a 3—level AMR calculation to compare with a uniform one for geodesic inflow to see the 
behavior of AMR in this problem. In Fig. [21 we plot AMR and uniform runs on top of each 
other for density vs. radial Schwarzschild coordinate. We see that while the AMR solution 
has reached a steady state at t = 151M, the uniform solution has not reached a steady state 
by the same time. 


D. Circular Motion of Test Particles 

We will now simulate the circular motion of a fluid with the numerical code. To do 
this we set up a circular flow with negligible pressure in the equatorial plane, in which 
angular velocity at each radial direction r is the Keplerian value, Ea. IBlOD . This is called a 
Keplerian. The last stable circular orbit for a particle moving around a Schwarzschild black 
hole is at r = 6 M (M is mass of black hole) Therefore the gas or particles will fall into the 
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black hole if their radial position is less than 6 M. When their radial position is bigger than 
6 M, they should rotate in a circular orbit. Here we simulate this problem and compare the 
numerical solution with the analytic expectations. This problem tests the code with sources 
in the 0 direction. 

In order to simulate this problem, we choose the computational domain to be 3 M < 
r < 20 M and 0 < 0 < 27T. The computational domain is filled with constant density 
and pressureless gas, rotating in circular orbits with the Keplerian velocity and zero radial 
velocity. In Fig. |21 the radial velocities of the gas vs. radial coordinate at different times are 
plotted. It is numerically observed that the gas inside the last stable orbit, r = 6 M, falls 
into the black hole while gas outside the last stable orbit follows circular motion with the 
Keplerian velocity as we expect analytically. In Fig. we plot the density of the fluid vs. 
radial coordinate at different times to see the behavior of the disk. It is clear from that gas 
falls into the black hole for r < 6 M. So the numerical results from our code are consistent 
with the analytic expectations. 

IV. NUMERICAL MODELING OF THE ACCRETION DISKS IN THE 
SCHWARZSCHILD COORDINATE 

In the present work, we do not intend to make a very concrete model for particular object, 
and it is assumed a rather simple initial configuration of gas. At t — 0 computational domain 
filled with some negligible values for density and pressure with zero radial and angular 
velocity. The injected value from companion star from outher boundary of computational 
domain are; = 1.0, Pi n = 10~ 3 /7, v r = 0.01 and v^ which is the Keplerian fluid velocity. 
Computational domain is 3 M < r < 100 M and 0 < 0 < 2i r. Where M is the mass of black 
hole. Following problems solved are 2D modeling on an equatorial plane. 

The choice of boundary conditions at the inner and the outher numerical boundary is 
rather important. The boundary conditions are treated as follows. The fictitious cell are 
placed just outside of a boundary and we prescribe physical variables in the fictitious cells. 
Numerical fluxes on the boundary wall are computed by solving a Riemann problems between 
states. The following boundary conditions are used. The value of physical variables in the 
fictitious cells are the same as those in the neighboring interior cells at all the times. Radial 
velocities direction can be changed depend on inner or outher boundary condition used. 
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A. 7 = 1.05 

The first results of the spiral shock wave are given on the accretion disk around the black 
hole using the Schwarzschild metric. In all simulations, the black hole is at the center of 
computational domain and it is represented as a white hole in the graphics. First, gas is 
injected from outlier boundary to accrete an accretion disk with spiral shock wave and then 
the injection is stopped to see the behavior of accretion disk during the evolutions. Finally 
the numerical simulation is stopped when solution reaches to steady state. 

The density of the accretion disk is illustrated in Fig© at t = 19643M. In early time of 
simulation accretion disk is not in steady state yet but two-armed spiral shock is already 
created and dynamical structure of the disk is not changed any more. The only thing which 
changes during the simulation after the certain evolution time is the amplitude of density. 
It is also showed in Fig© that the mass of accretion disk during the process of gas from 
outer boundary by the companion star is in the steady states between the evolution times 
10000M — 19643M. The two-armed spiral shock is created for the case 7 = 1.05 during the 
injecting gas. and also this spiral arms and accretion disk go to steady state. 

In order to watch the behavior of the spiral shock on an accretion disk while no matter 
gets into accretion disk, the injected gas is stooped. It is numerically observed that initial 
structure of two-armed spiral shock wave is slightly changed because the hydrodynamical 
forces given by the shock, which is created by injected gas, is gone. In order to balanced 
the forces on the accretion disk, the spiral shock waves are kept, which is more compact 
than the one in the beginning of this simulations, around the black hole. This newly formed 
two-armed spiral shocks are plotted at t = 47526 M in Fig© These shocks are in the steady 
state and two-arms spiral shock almost 180° apart to each other. In Fig© we depict ID 
cut at fixed radial coordinate, r = 20 . 12 M for density, radial velocity, orbital velocity and 
pressure at t = 47526 M. Two-armed spiral shocks are clearly observed. Angular momentum 
transform is also seen in that graphic because the orbital velocities of spiral arms are less 
than the orbital velocity of accretion disk. 
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B. 7 = 1.2 

In this case, the same problem is solved for accretion disk except 7 = 1.2 which is more 
expressible than the 7 = 1.05. In order to accrete an accretion disk, the same injection is 
made from outlier boundary with the case for 7 = 1.05 and numerical simulation is run 
until solution reach to steady state. Afterall the injection is stooped and the behavior of 
accretion disk, and spiral shock wave are numerically observed. 

We do injection from outer boundary continuously to create an accretion disk around 
the black hole for 7 = 1.2. In FigEl the two-armed spiral shock waves are created. The 
solution is almost in steady state in FigEl Eventhough the solution is not in steady state, 
the structure of accretion disk is not change during the evolutions. The result for 7 = 1.2 
comparable with Makita et al. (2000). But there are some differences around the black hole 
because our code fully relativistic and it can give detail structure when the matter close to 
black hole. Another differences in my code and Makita et al. (2000) is that we use inflow 
boundary condition that is allow to gas fall into black hole but they use freezing boundary. 
The using different boundary close to black hole also makes big difference in the structure 
of spiral shock waves. 

In order to understand the how spiral arms behave after companion star is removed, the 
injecting gas from outlier boundary is stooped. After the injecting is not allowed any more 
at t = 17502 M, FigEI] displays numerical result at t = 19475 M that two-armed spiral shock 
waves are still kept because of tidal forces between companion star and black hole. Mass of 
the accretion during the whole process is given at Fig ^2 It reaches a maximum mass, which 
is almost called steady-state point, whenever companion star is removed mass of accretion 
disk stars to decrease and goes to steady-state point. That point has also two-armed spiral 
shock waves. 


V. CONCLUSION 

The spiral shock waves on an accretion disk are an important mechanism to transform the 
angular momentum. In this simulation it is concluded that accretion disk is rotating with 
sub-Keplerian velocity, which is carrying positive angular momentum, and spiral shock waves 
moving much slower than the disk and they are carrying the negative angular momentum. 
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When the disk material hits to spiral wave, the angular momentum of accretion disk is 
transform to out of accretion disk. The shock waves around the black hole are a mechanism 
of transforming gravitational energy to radiation energy which is observed by the different 
x—ray observatory satellite, such as Chandra. 

Spiral arms in an accretion disk are one of the mechanism to emit x— ray. We have 
looked at the spiral structure, which is created at around the black hole when matter falls 
from companion star to primary star. Spiral structure in an accretion disk is created under 
certain condition. To create spiral arms adiabatic index must be less or equal than 1.2. 

In this simulations, three different simulations are done for 7 = 1.1,1.2 and 5/3. These 
simulations show that we do not have any spiral arms for 7 = 5/3 (the graphic for his case 
did not put here) but for the others. These results are consistent with results from Makita 
et al. (2000). We have also watched behavior of accretion disk after companion star stop 
injection. The two-armed spiral structure kept during the evolution. It concludes that the 
spiral shocks are formed by tidal forces, not by the inflow, of which claim was posed by 
Bisikalo et al. (1998b). 

From the point of view of dependence 7 , it is concluded that spiral arms are more tightly 
in smaller 7 cases than larger ones. Lower 7 means cooler disk with larger Mach number of 
the flow. Our results are comparable with those of Makita et al. (2000) who solve Newtonian 
hydrodynamical equation. Our and their results are also agree while adiabatic index, 7 , is 
bigger than 1.2, two armed spiral shock wave is not created which is observed for 7 = 5/3. 
Even the accretion disk may not be formed. 
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APPENDIX A: THE ANALYTIC SOLUTION OF GEODESICS FLOWS 

Pressureless gas, also called dust, falling onto a black hole in the radial direction, called 
geodesic flow. It can be also called free falling gas because there are no pressure forces 
opposing the inward motion of the gas. 

We use the fact that the elements of the accreting fluid fall along geodesics to get the 
analytic solution . In axisymmetric, steady-state flows the binding energy per baryon hU t 
is conserved. Hence for dust particles, h = 1 and the gravitational binding energy U t will 
remain constant. Since = —1, v r (r ) is now determined in terms of input an constant 

Ut and the known metric functions. Note that Ug = U<p = 0. 

First, we start with the geodesic equation for free falling dust: 


VrJJ = 0 . 
= {U% + Y^W)^ = 0 . 

Eq. ED can be rewritten using U a = dx a /dX, 


(Al) 


— + T^lPU p = 0. (A2) 

Substituting the index a — 0 = t in Eq. El gives 

fjTjt 

— + I%piruP = 0. (A3) 

Most of the Christoffel symbols at Eq. El are zero, except Y\ r = r* t = \g tt d r g t t = 
|(2 M/{r{r — 2 M))). The substituting the non-zero Christoffel symbols into Ea. (IA3li gives 
us 


dU l _ 2 M ut 

dr r(r — 2 M) 

where dU l /d\ = ( dU t /dr)(dr/d \). 

After doing some straightforward integration, Eq. El goes to 


(A4) 
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U 1 = _ 

To compute the radial component of the fluid velocity of geodesic gas, we need to know 
the radial four velocity of the gas. In order to compute that, we use the normalization of 
four velocities, which is 



U ^ = -1 


u'luU 1 + u r j„Y = -i. 


(A6) 


We substitute Eq. 4H3) into Eq. (Hi , to 


U r = 


2 M 


(AT) 


Using the relations between the four and three velocities, U r = Wv r , U l = Wj \J 1 — (2 M/r), 
we get 




(AS) 


where v r is the velocity which is observed by an observer outside the horizon. 

Now, we compute the density using the continuity equation from Eq. © which is 


dtiV^Wp) + di(av r D) = 0. (A9) 

Since we are looking for a steady state solution the time derivative of variables is zero and 
we have only the radial derivative in Eq. m- Then the density equation becomes 


di{av r D) = 0. 

After doing integration of Eq. cm we get 


(A10) 


av r D = d, 


(All) 
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where d is an integration constant. Now D can be computed from Eos. (IA8I) and (lAllj) . 


D 


d 


Finally, we compute the density p from Eos. (1A8I) and (1A12I) 


(A 12 ) 


1 d 

W r 2 (^)*(l - ’ 


where W is the Lorentz factor and given by 


(A13) 


W 



(A14) 


APPENDIX B: THE ANALYTIC REPRESENTATION OF CIRCULAR MOTION 

OF A TEST PARTICLE 


In this Appendix we compute the angular velocity and circular velocity of a particle on 
a circular orbit in the Schwarzschild spacetime, analytically. 

Since Schwarzschild geometry is time independent and spherically symmetric, the con¬ 
served quantities can be determined by the trajectory of particles. Because of spherical 
symmetry, motion is always defined in at a single plane and we can choose this plane to be 
the equatorial plane (6 = vr/2). Then 6 is constant in that plane for the motion of particles 
and the 9 derivatives vanish. The components of the momentum (Schutz et al. 1985) are 


P t = 9 U Pt = 


mE 


(1-T) 

r dr 
p = m — 
dr 

p e = 0 


P+ = 9 H V4> = 


mL 


(Bl) 


where m, E , r and L are the mass of particle, total energy, proper time and angular 
momentum, respectively. Here E = —p t /m and L = p^/m. 
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Now, we can derive the equation of motion for a particle in the equatorial plane using 
Eq. (IB2|) and the conservation relation, p.p = —m 2 . This gives us 



Eq. m can be rewritten by defining an effective potential V (r) and we get 

= E 2 - V 2 (r), (B3) 

where V' 2 W = ( 1-^)(1 + U). 

Ea. (IB3l) implies that since the left side of that equation is positive or zero, the total energy 
of a trajectory can be bigger or equal to the effective potential. 

In order to compute the angular velocity and period of a particle in a circular orbit, we 
differentiate Eq. (Ha with respect to r, and get 


d 2 r 

dr 2 


1 dV 2 {r) 

2 dr 


(B4) 


It is clear from Ea. (IB4l) that a circular orbit, which has constant r, is possible only at a 
minimum or maximum of the effective potential, V 2 {r). In a circular orbit r is constant and 
the left side of Ea. dBdl) goes to zero. If we take d 2 r/dr 2 = 0 and substitute in the expression 
for the effective potential, we can compute the circular orbit radius as 


L 2 
2 M 



(B5) 


From Ea. (IB5l) . a stable circular orbit at radius r has angular momentum which is 


L 2 = 


Mr 

2 _ 3 M ' 

r 


The total energy in a circular orbit is E~ — V 2 and it is 


E = 



(B6) 


(B7) 
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Now, the non zero components of the four velocity of a particle in the plane are 



(B8) 


and 


dt 

dr 



E 


(B9) 


r 


We find the angular velocity by dividing Eq. (1B8[) by Eq. ()B9[) : 


d(j) d(j)/dr 


dt dt/dr 



(BIO) 


this is called the Keplerian angular velocity. 

Finally, we can compute the circular velocity of a particle using the definition of four velocity 


relativistic hydrodynamical equations given in Section |TT1 , Eq. (IB Hi) and Eq. (IB9I) . We get 
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r(M) 


FIG. 1: The analytic solutions (red solid lines) with the numerical solutions (black circles) using 
256 zones in the radial direction for thermodynamical variables, p, v , and v r . 

















density 


Uniform and 3-levels AMR 



FIG. 2: The numerical solutions from a 3—level AMR run(64 zones, red circles) and a uniform grid 
run(160 zones, black straight line) for density vs. radial coordinate. 



r( Radial Velocity) 



FIG. 3: Radial velocity vs. r is plotted. The radial velocity of fluid in the disk stays zero for 
r > 6 M, during the evolution, r = 6 M is called the last stable circular orbit in a Keplerian disk. 
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FIG. 4: Density vs. r is plotted. Density of the fluid is plotted at different times using different 
colors. The matter falls into the black hole while r is less than 6 M. 
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FIG. 5: Plotting the density in the r — <fi plane with color for 7 = 1.05. It is taken at t = 19643M 
and it is in steady state. Two-armed spiral shock wave is created. 
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Time(M) = time*0.49 1e-5(M/Mo)(sec.) 


gamma = 1.05 



FIG. 6: Mass of accretion disk vs. time is plotted during the hole evolution. The injection is 
stopped at maximum mass. 
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FIG. 7: Zooming the interesting part of accretion disk at t = 47526 M to see two-armed shock wave 
clearly. 








FIG. 8: Plotting the density, radial velocity, orbital velocity and pressure for accretion disk at 
t = 47526 M in fixed r = 20.12 M. 
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FIG. 9: Plotting the density in the r — plane with color for 7 = 1.2. It is taken at t = 17502 M. 
It is in the steady state and two-armed spiral shock wave is already created and the structure of 
disk does not change. 
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FIG. 10: Plotting the density in the r — 4> plane with color for 7 = 1.2. It is taken at t = 19475M 
after injected is stopped atf = 17502M . It is seen that two-armed spiral shock wave is still kept 
because they are created tidal forces on the accretion disk. 
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time(M) = time * 0.49le-5 (M/M_o) (sec.) 

Mass of accretion disk vs. time. Injected gas is stopped at maximum mass 



FIG. 11: Mass of accretion disk vs. time is plotted during the hole evolution. The injection is 
stopped at maximum mass. 



